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Abstract. A practical introduction to stochastic modelling of reaction-diffusion processes is 
presented. No prior knowledge of stochastic simulations is assumed. The methods are explained 
using illustrative examples. The article starts with the classical Gillespie algorithm for the stochastic 
modelling of chemical reactions. Then stochastic algorithms for modelling molecular diffusion are 
given. Finally, basic stochastic reaction-diffusion methods are presented. The connections between 
stochastic simulations and deterministic models are explained and basic mathematical tools (e.g. 
chemical master equation) are presented. The article concludes with an overview of more advanced 
methods and problems. 
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1. Introduction. There are two fundamental approaches to the mathematical 
modelling of chemical reactions and diffusion: deterministic models which are based 
on differential equations; and stochastic simulations. Stochastic models provide a more 
detailed understanding of the reaction-diffusion processes. Such a description is often 
necessary for the modelling of biological systems where small molecular abundances 
of some chemical species make deterministic models inaccurate or even inapplicable. 
Stochastic models are also necessary when biologically observed phenomena depend on 
stochastic fluctuations (e.g. switching between two favourable states of the system). 

In this paper, we provide an accessible introduction for students to the stochastic 
modelling of the reaction-diffusion processes. We assume that students have a basic 
understanding of differential equations but we do not assume any prior knowledge of 
advanced probability theory or stochastic analysis. We explain stochastic simulation 
methods using illustrative examples. We also present basic theoretical tools which 
are used for analysis of stochastic methods. We start with a stochastic model of a 
single chemical reaction (degradation) in Section I2.fi introducing a basic stochastic 
simulation algorithm (SSA) and a mathematical equation suitable for its analysis (the 
so-called chemical master equation). Then we study systems of chemical reactions in 
the rest of Section [2J presenting the Gillespie SSA and some additional theoretical 
concepts. We introduce new theory whenever it provides more insights into the par- 
ticular example. We believe that such an example-based approach is more accessible 
for students than introducing the whole theory first. In Section [3] we study SSAs 
for modelling diffusion of molecules. We focus on models of diffusion which are later 
used for the stochastic modelling of reaction-diffusion processes. Such methods are 
presented in Section 01 We also introduce further theoretical concepts, including the 
reaction-diffusion master equation, the Smoluchowski equation and the Fokker-Planck 
equation. We conclude with Sections[5]and[6]where more advanced problems, methods 
and theory are discussed, giving references suitable for further reading. 

The stochastic methods and the corresponding theory are explained using several 
illustrative examples. We do not assume a prior knowledge of a particular computer 
language in this paper. A student might use any computer language to implement the 
examples from this paper. However, we believe that some students might benefit from 
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our computer codes which were used to compute the illustrative results in this paper. 
The computer codes (in Matlab or Fortran) can be downloaded from the website 
http://www.maths.ox.ac.uk/cmb/Education/ which is hosted by the Centre for 
Mathematical Biology in the Mathematical Institute, University of Oxford. 

2. Stochastic simulation of chemical reactions. The goal of this section is 
to introduce stochastic methods for the modelling of (spatially homogeneous) systems 
of chemical reactions. We present the Gillespie SSA, the chemical master equation and 
its consequences [HI HI]. We start with the simplest case possible, that of modelling a 
single chemical reaction, in Sect ion |2~T1 We then study two simple systems of chemical 
reactions in Sections 12.21 and 12.31 

2.1. Stochastic simulation of degradation. Let us consider the single chem- 
ical reaction 

A -^-> (2.1) 

where A is the chemical species of interest and k is the rate constant of the reaction. 
The symbol denotes chemical species which are of no further interest in what fol- 
lows. The rate constant k in (|2.1[) is defined so that k dt gives the probability that a 
randomly chosen molecule of chemical species A reacts (is degraded) during the time 
interval [t, t + dt) where t is time and dt an (infinitesimally) small time step. In par- 
ticular, the probability that exactly one reaction (|2.ip occurs during the infinitesimal 
time interval [t, t + dt) is equal to A{t)k dt where we denote the number of molecules 
of chemical species A at time t simply as A(t). This notational convention will be 
used throughout the paper. 

Let us assume that we have n molecules of A in the system at time t = 0, 
i.e. ^4(0) = no. Our first goal is to compute the number of molecules A(t) for 
times t > 0. To do that, we need a computer routine generating random numbers 
uniformly distributed in the interval (0, 1). Such a program is included in many modern 
programming languages (e.g. function rand in Matlab): It generates a number r 6 
(0, 1), so that the probability that r is in a subinterval (a, b) C (0, 1) is equal to b — a 
for any a, b £ (0, 1), a < b. 

The mathematical definition of the chemical reaction (12. 1|) can be directly used to 
design a "naive" SSA for simulating it. We choose a small time step At. We compute 
the number of molecules A(t) at times t = iAt, i = 1,2,3, ... , as follows. Starting 
with t = and A(0) — no, we perform two steps at time t: 

(al) Generate a random number r uniformly distributed in the interval (0, 1). 
(bl) If r < A(t)kAt, then put A(t + At) = A(t)-1; otherwise, A(t + At) = A(t). 
Then continue with step (al) for time t + At. 

Since r is a random number uniformly distributed in the interval (0, 1), the probability 
that r < A(t)k At is equal to A(t)k At. Consequently, step (bl) says that the proba- 
bility that the chemical reaction (|2.ip occurs in the time interval [t,t + At) is equal to 
A(t)k At. Thus step (bl) correctly implements the definition of (|2.1[) provided that 
At is small. The time evolution of A obtained by the "naive" SSA (al)-(bl) is given 
in Figure 0(a) for k = 0.1 sec^ 1 , A(0) = 20 and At = 0.005 sec. We repeated the 
stochastic simulation twice and we plotted two realizations of SSA (al)-(bl). We 
see in Figure l2~TT a) that two realizations of SSA (al)-(bl) give two different results. 
Each time we run the algorithm, we obtain different results. This is generally true 
for any SSA. Therefore, one might ask what useful and reproducible information can 
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Fig. 2.1. Stochastic simulation of chemical reaction \2. ijl /or k = 0.1 sec -1 and A(0) = 20. (a) 
Number of molecules of A as a function of time for two realizations of the "naive" SSA (al)-(bl) 
for At = 0.005 sec; (b) results of ten realizations of SSA ( a2)-( c2) (solid lines; different colours show 
different realizations) and stochastic mean \2.8t plotted by the dashed line. 

be obtained from stochastic simulations? This question will be addressed later in this 
section. 

The probability that exactly one reaction (|2.ip occurs during the infinitesimal time 
interval [t,t + dt) is equal to A(t)kdt. To design the SSA (al)-(bl), we replaced dt by 
the finite time step At. In order to get reasonably accurate results, we must ensure 
that A(t)kAt <C 1 during the simulation. We used k = 0.1 sec -1 and At = 0.005 sec. 
Since A(t) < A(0) = 20 for any t > 0, we have that A(t)k At e [0,0.01] for any 
t > 0. Consequently, the condition A(t)k At -C 1 is reasonably satisfied during the 
simulation. We might further increase the accuracy of the SSA (al)-(bl) by decreasing 
At. However, decreasing At increases the computational intensity of the algorithm. 
The probability that the reaction (|2.1[) occurs during the time interval [t, t + At) is 
less or equal to 1% for our parameter values. During most of the time steps, we 
generate a random number r in step (al) to find out that no reaction occurs in step 
(bl). Hence, we need to generate a lot of random numbers before the reaction takes 
place. Our next task will be to design a more efficient method for the simulation of 
the chemical reaction (|2.ip . We will need only one random number to decide when 
the next reaction occurs. Moreover, the method will be exact. There will be no 
approximation in the derivation of the following SSA (a2)-(c2). 

Suppose that there are A(t) molecules at time t in the system. Our goal is to com- 
pute time t + r when the next reaction (|2.1[) takes place. Let us denote by f(A(t), s) ds 
the probability that, given A(t) molecules at time t in the system, the next reaction 
occurs during the time interval [t + s, t + s + ds) where ds is an (infinitcsimally) small 
time step. Let g(A(t), s) be the probability that no reaction occurs in interval [t, t+s). 
Then the probability f(A(t), s) ds can be computed as a product of g(A(t),s) and the 
probability that a reaction occurs in the time interval [t + s, t + s + ds) which is given 
according to the definition of (|2.ip by A(t + s)kds. Thus we have 

f(A(t), s) ds = g{A{t),s)A{t + s)k ds. 

Since no reaction occurs in [t, t + s), we have A(t + s) = A(t). This implies 



f(A(t), s) ds = g(A(t),s)A(t)k ds. 



(2.2) 
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To compute the probability g(A(t), s), let us consider a > 0. The probability that no 
reaction occurs in the interval [t,t + a + da) can be computed as the product of the 
probability that no reaction occurs in the interval [t, t + a) and the probability that 
no reaction occurs in the interval [t + cr,t + a + da) . Hence 

g(A(t), a + da) = g(A(t), a)[l - A(t + a)k da}. 

Since no reaction occurs in the interval [t, t + a), we have A(t + a) = A(t). Conse- 
quently, 

9 (MtU + da)- 9 (A(t),a) = 
da 

Passing to the limit da — > 0, we obtain the ordinary differential equation (in the a 
variable) 

Solving this equation with initial condition g(A(t),0) = 1, we obtain 

g{A(t),a) = exp[-A(t)ka]. 

Consequently, (|2.2|) can be written as 

f(A(t),s)ds = A(t)kexp[-A(t)ks]ds. (2.3) 

Our goal is to find r such that t + r is the time when the next reaction occurs, 
provided that there are A(t) molecules of A in the system at time t. Such r 6 (0, oo) 
is a random number which has to be generated consistently with the definition of the 
chemical reaction (|2.1[) . To do that, we consider the function F(-) defined by 

F(t) = exp[-A(t)kT}. (2.4) 

The function F(-) is monotone decreasing for A{t) > 0. If r is a random number 
from the interval (0, oo), then F(t) is a random number from the interval (0,1). If 
r is a random number chosen consistently with the reaction (|2.1I) . then F(t) is a 
random number uniformly distributed in the interval (0, 1) which can be shown as 
follows. Let a, b, a < b, be chosen arbitrarily in the interval (0, 1). The probability 
that F(t) € (a, b) is equal to the probability that r S (F^ 1 ^), F^ 1 (a)) which is given 
by the integral of f(A(t),s) over s in the interval (F _1 (&), F^ 1 (a)). Using (|2.3p and 
(pT4j) . we obtain 

/ f(A(t),s)ds= / A{t)kexp[-A(t)ks]ds 

JF-l(b) JF- 1 ^) 



f F ~ W AT? 

-/ ^ds = -F[F- 1 (a)]+F[F- 1 {b)]=b-a. 

JF-Hb) ds 



Hence we have verifed that F(t) is a random number uniformly distributed in (0, 1). 
Such a number can be obtained using the random number generator (e.g. function 
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rand in Matlab). Let us denote it by r. The previous observation implies that we can 
generate the time step r by putting r = F(r). Using (|2.4p . we obtain 

r = exp[— A(t)kT]. 

Solving for r, we obtain the formula 

"1" 



In 



A(t)k*~ [r\ ■ (2 ' 5) 

Consequently, the SSA for the chemical reaction (|2.1|) can be written as follows. 
Starting with t — and A(0) — hq, we perform three steps at time t: 

(a2) Generate a random number r uniformly distributed in the interval (0, 1). 
(b2) Compute the time when the next reaction (|2.ip occurs as t + t where r is 
given by (|2.5|) . 

(c2) Compute the number of molecules at time t + r by A(t + r) = A(t) — 1. 
Then continue with step (a2) for time t + r. 

Steps (a2)-(c2) are repeated until we reach the time when there is no molecule of A 
in the system, i.e. A = 0. SSA (a2)-(c2) computes the time of the next reaction 
t + t using formula (12. 5|) in step (b2) with the help of one random number only. 
Then the reaction is performed in step (c2) by decreasing the number of molecules 
of chemical species A by 1. The time evolution of A obtained by SSA (a2)-(c2) is 
given in Figure I2~l7 b). We plot ten realizations of SSA (a2)-(c2) for k = 0.1 sec -1 
and A(0) — 20. Since the function A(t) has only integer values {0, 1,2,..., 20}, it is 
not surprising that some of the computed curves A(t) partially overlap. On the other 
hand, all ten realizations yield different functions A(t). Even if we made millions of 
stochastic realizations, it would be very unlikely (with probability zero) that there 
would be two realizations giving exactly the same results. Therefore, the details of 
one realization A(t) are of no special interest (they depend on the sequence of random 
numbers obtained by the random number generator). However, averaging values of 
A at time t over many realizations (e.g. computing the stochastic mean of A), we 
obtain a reproducible characteristic of the system - see the dashed line in Figure 
I2.1f b). The stochastic mean of A{t) over (infinitely) many realizations can be also 
computed theoretically as follows. 

Let us denote by p n (t) the probability that there are n molecules of A at time t 
in the system, i.e. A(t) — n. Let us consider an (infinitesimally) small time step dt 
chosen such that the probability that two molecules are degraded during [t, t + dt) 
is negligible compared to the probability that only one molecule is degraded during 
[t, t + dt). Then there are two possible ways for A(t + dt) to take the value n: either 
A(t) = n and no reaction occurred in [t, t + dt), or A(t) = n + 1 and one molecule was 
degraded in [t, t + dt). Hence 

p n (t + dt) = p n (t) x (1 — kndt) + p n+ i(t) x k(n + 1) dt. 

A simple algebraic manipulation yields 

p n (t + dt) -p n (t) 



dt 



k(n + l)p n+1 (t) - knp n {t). 



Passing to the limit dt — > 0, we obtain the so-called chemical master equation in the 
form 

= k(n + l)p n+1 - knp n . (2.6) 

dt 
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Equation (12. 6p looks like an infinite system of ordinary differential equations (ODEs) 
for p n , n = 0, 1, 2, 3, ... . Our initial condition A(0) = n implies that there are never 
more than no molecules in the system. Consequently, p n = for n > no and the 
system (|2.6j) reduces to a system of (no + 1) ODEs for p n , n < n . The equation for 
n = no reads as follows 



dt 



-kn p no . 



Solving this equation with initial condition p no (0) = 1, we get p no (t) = exp[— fcriot]. 
Using this formula in the chemical master equation (|2.6[) for p„ _i(t), we obtain 

j£Pno-l(*) = kn a exp[-kn t] - k(n a - l)p„ _i(t). 

Solving this equation with initial condition p„ o _i(0) = 0, we obtain p„ _i(i) = 
exp[— k(no — l)i]no(l — exp[— kt]). Using mathematical induction, it is possible to 
show 

Pn (t) = cxp[-knt] M {1 - exp[-kt]} na ~ n . (2.7) 

The formula (|2.7p provides all the information about the stochastic process which 
is defined by (|2.1[) and initial condition A(0) = no. We can never say for sure that 
A(t) — n; we can only say that A(t) = n with probability p n (t). In particular, formula 
(|2.7p can be used to derive a formula for the mean value of A(t) over (infinitely) many 
realizations, which is defined by 



M{t) = J^n Pn (t). 

n=Q 

Using (|2.7p . we deduce 

n n Q , s 

M(t) = J2 n Pn(t) = J2 nexp[-knt] I "° Ul- exp[-fcf]} n °"" 

/n 0-A fl r (no-l)-(n-l) r r L.n™" 1 



n=0 



= n exp[-fet]^ r°_ M{l-exp[-iW]} lBo " 1, - lT, - i; {ex P [-jW]} T 
= n exp[-fci]. (2.8) 



The chemical master equation (|2.6p and its solution (|2.7p can be also used to quan- 
tify the stochastic fluctuations around the mean value (|2.8|) . i.e. how much can an 
individual realization of SSA (a2)-(c2) differ from the mean value given by (|2.8p . We 
will present the corresponding theory and results on a more complicated illustrative 
example in Section l2.2l Finally, let us note that a classical deterministic description of 
the chemical reaction (|2.ip is given by the ODE da/dt — —ka. Solving this equation 
with initial condition a(0) = no, we obtain the function (|2.8p . i.e. the stochastic mean 
is equal to the solution of the corresponding deterministic ODE. However, we should 
emphasize that this is not true for general systems of chemical reactions, as we will 
see in Section 12.31 and Section 15.11 
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2.2. Stochastic simulation of production and degradation. We consider 
a system of two chemical reactions 

A 0, A. (2.9) 

The first reaction describes the degradation of chemical A with the rate constant k\. 
It was already studied previously as reaction (|2.ip . We couple it with the second 
reaction which represents the production of chemical A with the rate constant k 2 . 
The exact meaning of the second chemical reaction in (|2.9[) is that one molecule of 
A is created during the time interval [t, t + dt) with probability k 2 dt. As before, the 
symbol denotes chemical species which are of no special interest to the modeller. 
The impact of other chemical species on the rate of production of A is assumed to be 
time independent and is already incorporated in the rate constant k 2 - To simulate 
the system of chemical reactions (|2.9p . we perform the following four steps at time t 
(starting with A(0) = uq at time t = 0): 

(a3) Generate two random numbers r±, r 2 uniformly distributed in (0, 1). 
(b3) Compute a = A{t)k\ + k 2 . 

(c3) Compute the time when the next chemical reaction takes place as t + r where 



iln 

a 



(2.10) 

(d3) Compute the number of molecules at time t + r by 

A(t + T) = \ A V + ] i ! r2 J^ a ° ; (2.11) 
[A(t) — 1 it r 2 >k 2 /a - 

Then continue with step (a3) for time t + r. 

To justify that SSA (a3)-(d3) correctly simulates (|2.9p . let us note that the probability 
that any of the reactions in (|2.9p takes place in the time interval [t, t + dt) is equal to 
ao dt. It is given as a sum of the probability that the first reaction occurs (A(t)kidt) 
and the probability that the second reaction occurs (k 2 dt) . Formula (|2.10p gives the 
time t + t when the next reaction takes place. It can be justified using the same 
arguments as for formula (|2.5p . Once we know the time t + r, a molecule is produced 
with probability k 2 /ctQ, i.e. the second reaction in (|2.9p takes place with probability 
k 2 /ao- Otherwise, a molecule is degraded, i.e. the first reaction in (|2.9p occurs. The 
decision as to which reaction takes place is given in step (d3) with the help of the 
second uniformly distributed random number r 2 . 

Five realizations of SSA (a3)-(d3) are presented in Figurc [2~2T a) as solid lines. We 
plot the number of molecules of A as a function of time for A(0) = 0, k\ = 0.1 sec -1 
and k 2 = 1 sec -1 . We see that, after an initial transient, the number of molecules A(t) 
fluctuates around its mean value. To compute the stochastic mean and quantify the 
stochastic fluctuations, we use the chemical master equation which can be written for 
the chemical system (|2.9p in the following form 

— — = ki{n + l)p n +i ~ hnp n + k 2 p n -i - k 2 p n (2-12) 
dt 

where p n (t) denotes the probability that A(t) = n for n = 0, 1, 2, 3, Let us note 

that the third term on the right hand side is missing in (|2.12p for n = 0; i.e. we use 
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Fig. 2.2. Stochastic simulation of the system of chemical reactions \2.9) for A(0) = 0, k\ = 
0.1 sec -1 and ki = 1 sec - 1 . (a) A(t) given by five realizations of SSA (a3)-(d3) (solid lines) and 
stochastic mean (dashed line), (b) Stationary distribution <p( n ) obtained by long time simulation of 
SSA (a.3)-(d3) (gray histogram) and by formulae \2.21\ - \2.22\ (red solid line). 



the convention that p_i = 0. The first two terms on the right hand side correspond 
to the first reaction in (|2.9[) . They already appeared in equation (|2.6[) . Production 
of A is described by the third and fourth term on the right hand side of (|2.9|) . To 
derive the chemical master equation (|2.12p . we can use similar arguments as in the 
derivation of (|2.6p . The stochastic mean M(t) and variance V(t) are defined by 

oo oo 

M(t) =X>P»W, V(t) = Y,(n-M(t)) 2 Pn (t). (2.13) 

n=0 n=0 

The stochastic mean M(t) gives the average number of molecules of A at time t, 
while the variance V(t) describes the fluctuations. In Section 12. li we first solved 
the chemical master equation (|2.6|) and then we used its solution (|2.7|) to compute 
M(t). Alternatively, we could use the chemical master equation to derive an evolution 
equation for M(t), i.e. we could find M(t) without solving the chemical master 
equation. Such an approach will be presented in this section. Multiplying (|2.12[) by 
n and summing over n, we obtain 

j oo oo oo oo oo 

^2 n P n = ki X! n ( n + ~ fci X! n2pn + k2 *^2 n Pn~i - h y^^np n . 

n— n— n— n— 1 n— 

Using definition (|2.13p on the left hand side and changing indices n + 1 — > n (resp. 
n — 1 — > n) in the first (resp. third) sum on the right hand side, we obtain 

— - = ki 2J(«- l)np n - ki ^2n 2 p n + k 2 >^(n + l)p„ - A; 2 «Pn- 

n— n— n— n— 

Adding the first and the second sum (resp. the third and the fourth sum) on the right 
hand side, we get 

dM ,™ 7 ~ , n „. 

—^j- = -kl n Pn + k 2 2_^Pn- (2-14) 

ra=0 n=0 
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Since p n (t) is the probability that A(t) = n and A(t) is equal to a nonnegative integer 
with probability 1, we have 



X>n(*) = l- (2-15) 



n=0 

Using this fact together with the definition of M(t), equation (|2.14p implies the evo- 
lution equation for M(t) in the form 

™£ = _ klM + k2 _ (2.16) 
dr 

The solution of (|2.16|) with initial condition M(Q) = is plotted as a dashed line 
in Figure [2~2T a). To derive the evolution equation for the variance V(t), let us first 
observe that definition (|2.13[) implies 

oo 

Y / n 2 Pn{t) =V(t) + M(tf. (2.17) 

n=0 

Multiplying (j2~T2"]) by n 2 and summing over n, we obtain 

, oo oo oc oo oo 

-g, ^2 n 2 p n = ki ^ n2 ( n + l )Pn+i ~ k i^2 n3 Pn + k 2^2 ™ 2 Pn-i - k 2 X! n2 Pn- 

n— n— n— n— 1 n— 

Changing indices n + 1 — ► n (resp. n — 1 — ► n) in the first (resp. third) sum on the 
right hand side and adding the first and the second sum (resp. the third and the 
fourth sum) on the right hand side, we get 

^ oc oc oo 

— ^ n2 Pn = k l X^ 2 " 2 + U ^ Pn + fc2 ^2( 2n + l ^ Pn - 

n—0 71=0 71—0 

Using (|2Tf|) . (j2~T5| and (prig]) , we obtain 

— + 2M— = -2ki \V + M 2 ] + k x M + 2k 2 M + k 2 . 
at at 1 J 

Substituting (|2.16p for dM/dt, we derive the evolution equation for the variance V(t) 
in the following form 

dV 

— = -2kiV + kiM + k 2 . (2.18) 
at 

The time evolution of M{t) and V{t) is described by (j2~T6|) and (|2T8]l . Let us define 
the stationary values of M{t) and V{t) by 

M s = lim M(i), V s = lim V(t). (2.19) 

t — >oo t — >oo 

The values of M s and V s can be computed using the steady state equations corre- 
sponding to (|2.16p and (|2.18p . namely by solving 

Q = -kiM s + k 2 , and = -2k x V s + kiM s + fc 2 . 
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Consequently, 



M s = V s = k l. 

ki 



For our parameter values fci = 0.1 sec -1 and fc 2 = lsec -1 , we obtain M s = V s = 10. 
We see in Figure [2~2T a) that A{t) fluctuates after a sufficiently long time around the 
mean value M s = 10. To quantify the fluctuations, one often uses the square root of 
V s , the so-called mean standard deviation which is equal to \/l0. 

More detailed information about the fluctuations is given by the so-called sta- 
tionary distribution 4>{n), n = 0, 1, 2, 3, . . . , which is defined as 

Mri) = lim pJt). (2.20) 

t — >oo 

This means that 4>{n) is the probability that A(t) = n after an (infinitely) long 
time. One way to compute 4>{n) is to run SSA (a3)-(d3) for a long time and create 
a histogram of values of A(t) at given time intervals. Using fci = 0.1 sec -1 and 
fc 2 = lscc -1 , the results of such a long time computation are presented in Figure 
12.2( b) as a gray histogram. To compute it, we ran SSA (a3)-(d3) for 10 5 seconds, 
recording the value of A(t) every second and then dividing the whole histogram by 
the number of recordings, i.e. by 10 5 . An alternative way to compute 4>(n) is to use 
the steady state version of the chemical master equation (12.12[) . namely 

= fci 0(1) - fc 2 0(0) 

= fci (ri + 1) 0(ri + 1) — kin <j){n) + fc 2 (j)(n — 1) — fc 2 </>(ro), f° r n ^ 1> 
which implies 

4>{n + l) = — : -\kw<f>(ri) + fc 2 0(n) - fc 2 0(n- 1)1, forn>l. (2.22) 

fei(n + 1) 

Consequently, we can express cf>(n) for any n > 1 in terms of </>(0). The formulae 
(|2.2ip - (|2.22p yield an alternative way to compute <p(n). We put 0(0) = 1 and we 
compute 0(n), for sufficiently many n, by (|2.2ip - (|2.22p . Then we divide <fi(n), n > 0, 
by S0(n). The results obtained by ([2"^ - (f2~2"2"]) are plotted in Figure ^T[b) as a 
(red) solid line. As expected, the results compare well with the results obtained by 
the long time stochastic simulation. 

We can also find the formula for <fi(n) directly. We let a reader to verify that the 
solution of the recurrence formula (|2.2ip - (|2.22p can be written as 

*(») = §(irf ( 2 - 23 ) 



n! \ fci 

where C is a real constant. Using (|2.15[) and (|2.20p . we have 



£>(») = !. (2.24) 



n=0 



Substituting (|2.23|1 into the normalization condition (|2.24[) . we get 



(k 2 y r ~ i /fc 2 y 



n\ V fci J ^— i n\ V fci 

n=0 x 7 n=0 v 
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where we used the Taylor series for the exponential function to get the last equal- 
ity. Consequently, C = exp[— k 2 /k{\ which, together with (12.23p . implies that the 
stationary distribution <j)(n) is the Poisson distribution 



1 

n! 





' k 2 


(1) exp 





(2.25) 



Thus the red solid line in Figure I2.2f b) which was obtained numerically by the re- 
currence formula (|2.21|) — (|2.22|) can be also viewed as the stationary distribution <j>(n) 
given by the explicit exact formula (|2.25p . 

2.3. Gillespie algorithm. SSAs (a2)-(c2) and (a3)-(d3) were special forms of 
the so-called Gillespie SSA. In this section, we present this algorithm for a more com- 
plicated illustrative example which will also involve second-order chemical reactions. 
Such chemical reactions are of the following form 



A + A^C, A + B^D. (2.26) 

In the first equation, two molecules of A react with rate constant k\ to produce 

C . The probability that the reaction takes place in the time interval [t, t + dt) is 
equal to A(t)(A(t) — l)fcidt. We define the propensity function of the first reaction as 
ai(t) = A(t)(A(t) — l)k\. Then the probability that the first reaction occurs in the 
time interval [t, t + dt) is equal to a±(t) dt. The propensity function which corresponds 
to the second equation in (|2.26|) is defined as ct2(t) = A{t)B(f)k\ and the probability 
that the second reaction occurs in the time interval [t, t + dt) is equal to ct2(t) dt. In 
such a case, one molecule of A and one molecule of B react to form a molecule of 

D. In general, the propensity function can be defined for any chemical reaction so 
that its product with dt gives the probability that the given reaction occurs in the 
innnitesimally small time interval [t, t + dt). 

We consider that A and B can react according to l|2.26[) . Moreover, we assume 
that they are also produced with constant rates, that is, we consider a system of four 
chemical equations: 

A + A^>® A + B^V) (2.27) 



A B. (2.28) 

Let us note that we are not interested in chemical species C and D. Hence, we 
replaced them by 0, consistent with our previous notation of unimportant chemical 
species. To simulate the system of chemical reactions (|2.27p - (|2.28[) . we perform the 
following four steps at time t (starting with A(0) — no, -B(O) = mo at time t = 0): 

(a4) Generate two random numbers r%, r 2 uniformly distributed in (0, 1). 
(b4) Compute the propensity functions of each reaction by a-y = A(t)(A(t) — l)fei, 
ct2 — A{t)B{t)k 2 , c*3 = &3 and = foj,. Compute ao = a± + a 2 + 0:3 + 014. 
(c4) Compute the time when the next chemical reaction takes place as t+T where 



1 

n 



(2.29) 
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(a) 



(b) 




40 60 
time [sec] 




40 60 
time [sec] 



Fig. 2.3. Five realizations of SSA (a4)-(d4). Number of molecules of chemical species A 
(left panel) and B (right panel) are plotted as functions of time as solid lines. Different colours 
correspond to different realizations. The solution of i2. 331 34\ is given by the dashed line. We 
use A(0) = 0, B(0) = 0, fci = 10 -3 sec' 1 , k 2 = 10 -2 sec' 1 , k 3 = 1.2 sec' 1 and fc 4 = 1 sec' 1 . 



(d4) Compute the number of molecules at time t + r by 

{A(t)-2 if < r 2 < ai/ao; 

A(t)-\ if a-i/ao < r 2 < («i + a 2 )/a ; 

A(t) + 1 if (ai + a 2 )/a Q < r 2 < (ai + a 2 + a 3 )/a Q ; 
A(t) if (ai + a 2 + a 3 )/aa < r 2 < 1; 



(2.30) 



B(t + r) = 



B(t) 
B(t) - 1 

B{t) 
B(t) + 1 



if < r 2 < ai/a ; 
if ai/a < r 2 < [a\ + a 2 )/a ; 
if (ai + a 2 )/a < r 2 < {a\ + a 2 + a 3 )/a ; 
if [a\ + a 2 + a 3 )/a < r% < 1; 



(2.31) 



Then continue with step (a4) for time t + r. 



SSA (a4)-(d4) is a direct generalisation of SSA (a3)-(d3). At each time step, we 
first ask the question when will the next reaction occur? The answer is given by 
formula (|2.29p which can be justified using the same arguments as formulae (|2.5[) or 
(|2.10p . Then we ask the question which reaction takes place. The probability that 
the i-th chemical reaction occurs is given by on/a®. The decision which reaction takes 
place is given in step (d4) with the help of the second uniformly distributed random 
number r 2 . Knowing that the i-th reaction took place, we update the number of 
reactants and products accordingly. Results of five realizations of SSA (a4)-(d4) are 
plotted in Figure [23] as solid lines. We use A(0) = 0, B(0) = 0, fci = lO^^ec" 1 , 
k 2 = 10~ 2 sec -1 , fc 3 = 1.2 sec -1 and £4 = lsec -1 . We plot the number of molecules 
of chemical species A and B as functions of time. We see that, after initial transients, 
A(t) and B(t) fluctuate around their average values. They can be estimated from 
long time stochastic simulations as 9.6 for A and 12.2 for B. 

Let p n ,m(t) be the probability that A(t) = n and B(t) — m. The chemical master 
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equation can be written in the following form 

^?' m = h(n + 2)(n + l)p n+2 ,m - hn(n - l)p n , m 
+ k 2 (n + l)(m + l)p n +i,m+i - k 2 nmp n , m 

+ k 3 Pn-l,7n ~ k 3 p n ^ m + /c 4 p„ !m _i - /C4 Pn,m (2.32) 

for 7i, ?7i ^ 0, with, the convention that Pn,m — 

if n < or m < 0. The first 
difference between (|2.32p and the chemical master equations from the previous sections 
is that equation (|2.32p is parametrised by two indices n and m. The second important 
difference is that (|2.32p cannot be solved analytically as we did with (|2.6[) . Moreover, 
it does not lead to closed evolution equations for stochastic means and variances; i.e. 
we cannot follow the same technique as in the case of equation (|2 . 1 2[) . The approach 
from the previous section does not work. Let us note that the probability p n ,m(t) is 
sometimes denoted by p(n,m,t); such a notational convention is often used when we 
consider systems of many chemical species. We will use it in the following sections to 
avoid long subscripts. 

The classical deterministic description of the chemical system (|2.27[) - (12.28p is 
given by the system of ODEs 

^ = -2fcia 2 - k 2 ab + k 3 , (2.33) 

— = -k 2 ab + k i . (2.34) 
at 

The solution of (|2~35|l - (|2~34"|l with initial conditions o(0) = and 6(0) = is plotted 
as a dashed line in Figure [231 Let us note that the equations (|2.33p ~ (|2.34p do not 
describe the stochastic means of A(t) and B(t). For example, the steady state values 
of (|2.33p - (|2.34p are (for the parameter values of Figure l2.3p equal to a s = b s = 10. On 
the other hand, the average values estimated from long time stochastic simulations are 
9.6 for A and 12.2 for B. We will see later in Section l5~Tl that the difference between 
the results of stochastic simulations and the corresponding ODEs can be even more 
significant. 

The stationary distribution is defined by 

4>{n,m) = lim p n>m (t). 

This can be computed by long time simulations of SSA (a4)-(d4) and is plotted in 
Figure [2T4ja). We see that there is a correlation between the values of A and B. This 
can also be observed in Figure [2T751 Looking at the blue realizations, we see that the 
values of A(t) are below the average and the values of B{t) are above the average, 
similarly for other realizations. One can also define the stationary distribution of A 
only by 



5^(n,m). (2.35) 



m=0 



Summing the results of Figure l2.4r a) over m, we obtain <p(n) which is plotted in Figure 
2.4f b) as a gray histogram. The red bar highlights the steady state value a s = 10 of 
system (|2T33]l - (f2T34|l . 

SSAs (a3)-(d3) and (a4)-(d4) were special forms of the so-called Gillespie SSA. 
To conclude this section, we formulate the Gillespie SSA in its full generality. Let us 
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(a) (b) 




number of A molecules number of A molecules 



Fig. 2.4. (a) Stationary distribution <j>(n,m) obtained by long time simulation of (a4)-(d4) for 
ki = 10 -3 sec" 1 , &2 = 10 — 2 sec -1 , kg = 1.2 sec -1 and kn = 1 sec -1 , (b) Stationary distribution of 
A obtained by i2.35l . 



consider that we have a system of q chemical reactions. Let ai(t) be the propensity 
function of the i-th reaction, i = 1, 2, . . . , q, at time t, that is, a^(i) dt is the probability 
that the i-th reaction occurs in the time interval [t,t + dt). Then the Gillespie SSA 
consists of the following four steps at time t. 

(a5) Generate two random numbers r±, r 2 uniformly distributed in (0, 1). 
(b5) Compute the propensity function on(t) of each reaction. Compute 

1=1 

(c5) Compute the time when the next chemical reaction takes place as t + r where 

t is given by (|2~^|) . 
(d5) Compute which reaction occurs at time t + t. Find j such that 

1 • 7 ~ 1 1 J 

T2 > — 7 a, and r 2 < — > «i . 

i—i i=i 

Then the j-th reaction takes place, i.e. update numbers of reactants and 
products of the j-th reaction. 
Continue with step (a5) for time t + r. 

The Gillespie SSA (a5)-(d5) provides an exact method for the stochastic simulation 
of systems of chemical reactions. It was applied previously as SSA (a2)-(c2) for 
the chemical reaction (|2.1|) . as SSA (a3)-(d3) for the chemical system (|2.9|) and as 
SSA (a4)-(d4) for the chemical system (|2.27p ~ (|2.28p . Our simple examples can be 
simulated quickly in Matlab (in less than a second on present-day computers). If 
one considers systems of many chemical reactions and many chemical species, then 
SSA (a5)-(d5) might be computationally intensive. In such a case, there are ways 
to make the Gillespie SSA more efficient. For example, it would be a waste of time 
to recompute all the propensity functions at each time step (step (b5)). We simu- 
late one reaction per one time step. Therefore, it makes sense to update only those 




propensity functions which arc changed by the chemical reaction which was selected in 
step (d5) of SSA (a5)-(d5). A more detailed discussion about the efficient computer 
implementation of the Gillespie SSA can be found e.g. in [16J. 

3. Diffusion. Diffusion is the random migration of molecules (or small particles) 
arising from motion due to thermal energy [3] . As shown by Einstein, the kinetic en- 
ergy of a molecule (e.g. protein) is proportional to the absolute temperature. In 
particular, the protein molecule has a non-zero instantaneous speed at, for exam- 
ple, room temperature or at the temperature of the human body. A typical protein 
molecule is immersed in the aqueous medium of a living cell. Consequently, it can- 
not travel too far before it bumps into other molecules (e.g. water molecules) in the 
solution. As a result, the trajectory of the molecule is not straight but it executes 
a random walk as shown in Figure l3~TT a). We plot six possible trajectories of the 
protein molecule with six different colours. All trajectories start at the origin and are 
followed for 10 minutes. We will provide more details about this figure together with 
the methods for simulating molecular diffusion in the rest of this section. Stochastic 
models of diffusion which are based on the Smoluchowski equation are introduced in 
Section 13.11 In Section 13. 2\ we introduce a model which is suitable for coupling with 
the Gillespie SSA. Both modelling approaches will be used later in Section [1] for the 
stochastic modelling of reaction-diffusion processes. Let us note that there exist other 
models of molecular diffusion - they will be discussed in Section [6l 

3.1. Smoluchowski equation and diffusion. Let [X(t), Y(t), Z(t)] € K 3 be 
the position of a diffusing molecule at time t. Starting with [X(0), Y(0), Z(0)] = 
[xq, yo, zq], we want to compute the time evolution of [X(t), Y(t), Z(t)]. To do that, we 
make use of a generator of random numbers which are normally distributed with zero 
mean and unit variance. Such a generator is part of many modern computer languages 
(e.g. function randn in Matlab). Diffusive spreading of molecules is characterised by 
a single diffusion constant D which depends on the size of the molecule, absolute 
temperature and viscosity of the solution [3]. Choosing time step At, we compute 
the time evolution of the position of the diffusing molecule by performing two steps 
at time t: 
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(a6) Generate three normally distributed (with zero mean and unit variance) 

random numbers £ x , £ y and £ z . 
(b6) Compute the position of the molecule at time t + At by 

X(t + At) =X(t)+ V2D At £ x , (3.1) 
Y(t + At) = Y(t) + \/2D At £ y , (3.2) 
Z(t + At) = Z{t) + V2D At (3.3) 

Then continue with step (a6) for time t + At. 

Choosing D = 1CP 4 mm 2 sec -1 (diffusion constant of a typical protein molecule), 
[X(0),Y(0),Z(0)] = [0,0,0] and At = 0.1 sec, we plot six realizations of SSA (a6)- 
(b6) in Figure l3~TT a). We plot only the x and y coordinates. We follow the diffusing 
molecule for 10 minutes. The position of the molecule at time t = 10 min is denoted 
as a black circle for each trajectory. 

Equations (|3.ip - (|3.3[) are discretized versions of stochastic differential equations 
(SDEs) which are sometimes called Smoluchowski equations. Some basic facts about 
SDEs can be found e.g. in [21 HI]- A more accessible introduction to SDEs can be 
found in 23J which has a similar philosophy as our paper. Reference [23] is a nice 
algorithmic introduction to SDEs for students who do not have a prior knowledge 
of advanced probability theory or stochastic analysis. We will not go into details of 
SDEs in this paper, but only highlight some interesting facts which will be useful 
later. 

First, equations (|3. 1[> — (|3.3[) are not coupled. To compute the time evolution of 
X(t), we do not need to know the time evolution of Y(t) or Z(t). We will later focus 
only on the time evolution of the x-th coordinate, effectively studying one-dimensional 
problems. Two-dimensional or three-dimensional problems can be treated similarly. 
Second, we see that different realizations of SSA (a6)-(b6) give different results. To get 
more reproducible quantities, we will shortly study the behaviour of several molecules. 
However, even in the case of a single diffusing molecule, there are still quantities 
whose evolution is deterministic. Let <p(x, y, t) dxdydz be the probability that X(t) £ 
[x, x + dx), Y(t) £ [y, y + dy) and Z(t) £ [z,z + dz) at time t. It can be shown that 
if evolves according to the partial differential equation 

dip _ { d 2 ip d 2 ip t d 2 ip 
dt \ dx 2 dy 



+ ^ + m (3-4) 



which is a special form of the so-called Fokker-Planck equation. Since our random walk 
starts at the origin, we can solve ([3.4| with initial condition tp(x, y, z,0) = 5(x,y,z) 
where S is the Dirac distribution at the origin. We obtain 



f(x,y,z,t) = (4p7rf)3/2 exp 



y 2 + z 2 



4Dt 



In order to visualise this probability distribution function, we integrate it over z to 
get probability distribution function 



( •(.>■. ii. I) : : / <p(x, y, z, t)dz = — - exp 



x 2 + y 2 
ADt 



(3.5) 



This means that ip(x, y, t) dxdy is the probability that X(t) £ [x, x + dx) and Y(t) £ 
[y,y + dy) at time t. The function ij}{x,y,t) at time t = 10 min is plotted in Figure 
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I3.1f b). It can be obtained also by computing many realizations of SSA (a6)-(b6) and 
plotting the histogram of positions of a molecule at time 10 min; such positions were 
denoted as black circles for the six illustrative trajectories in Figure [3~TT a). 

One important issue which was not addressed previously is that molecules diffuse 
in bounded volumes, i.e. the domain of interest has boundaries and suitable bound- 
ary conditions must be implemented. In the rest of this paper, we focus on one- 
dimensional problems to avoid technicalities. Hence, we effectively study diffusion of 
molecules in the one-dimensional interval [0,L]. Then the SSA can be formulated as 
follows: 

(a7) Generate a normally distributed (with zero mean and unit variance) random 
number £. 

(b7) Compute the position of the molecule at time t + At by 

X(t + At) = X(t) + V2D At f. (3.6) 

(c7) If X(t + At) computed by (|3~6)) is less than 0, then 
X(t + At) = -X(t) - V2D At f. 

If X(t + At) computed by (|3.6| is greater than L, then 

X(t + At) = 2L- X{t) - V2DAt £. 

Then continue with step (a7) for time t + At. 

The boundary condition implemented in step (c7) is the so-called reflective boundary 
condition or zero flux boundary condition. It can be used when there is no chem- 
ical interaction between the boundary and diffusing molecules. More complicated 
boundary conditions are discussed in [3 [8] . 

Choosing D = 1CP 4 mm 2 sec" 1 , L = 1 mm, X(0) = 0.4 mm and At = 0.1 sec, we 
plot ten realizations of SSA (a7)-(c7) in Figure I3.2f a) . Let us assume that we have 
a system of 1000 molecules which are released at position x — 0.4 mm at time t = 0. 
Then Figure 13.21( a) can be viewed as a plot of the trajectories of ten representative 
molecules. Considering 1000 molecules, the trajectories of individual molecules are 
of no special interest. We are rather interested in spatial histograms (density of 
molecules). An example of such a plot is given in Figure l3"T2T b'). We simulate 1000 



18 



RADEK ERBAN ET AL. 



molecules, each following SSA (a7)-(c7). At time t = 4 min, we divided the domain 
of interest [0, L\ into 40 bins of length h = L/40 = 25 /j,m. We calculated the number 
of molecules in each bin [(i — l)h, ih), i = 1, 2, . . . , 40, at time t = 4 min and plotted 
them as a histogram. 

Let us note that the deterministic counterpart to the stochastic simulation is a 
solution of the corresponding Fokker-Planck equation (diffusion equation in our case) 
which, in one dimension with zero flux boundary conditions, reads as follows 

The solution of (|3.7p with the Dirac-like initial condition at x = 0.4 mm is plotted as 
a red solid line in Figure 13.2( b) for comparison. 



3.2. Compartment-based approach to diffusion. In Section 13. 1[ we sim- 
ulated the behaviour of 1000 molecules by computing the individual trajectories of 
every molecule (using SSA (a7)-(c7)). At the end of the simulation, we divided the 
computational domain [0, L] into K = 40 compartments and we plotted numbers of 
molecules in each compartment in Figure [3T2^b) . In particular, most of the com- 
puted information (1000 trajectories) was not used for the final result - the spatial 
histogram. We visualised only 40 numbers (numbers of molecules in compartments) 
instead of 1000 computed positions of molecules. In this section, we present a different 
SSA for the simulation of molecular diffusion. We redo the example from Section |3~T1 
but instead of simulating 1000 positions of the individual molecules, we are going to 
simulate directly the time evolution of 40 compartments. 

To do that, we divide the computational domain [0, L] into K = 40 compartments 
of length h = L/K = 25 /an. We denote the number of molecules of chemical species 
A in the z-th compartment [(i — l)h, ih) by A*, i = 1, . . . , K . We apply the Gillespie 
SSA to the following chain of "chemical reactions" : 

d d d d 

A 2 ^ A 3 ^ ... i=i A K (3.8) 

d d d d 

where 

d 

Ai < — A; + i means that Ai — > Aj+x and Aj+i — > Ai. 

d 

We will shortly show that the Gillespie SSA of (13. 8|) provides a correct model of 
diffusion provided that the rate constant d in (|3.8p is chosen as d = D/h 2 where 
D is the diffusion constant and h is the compartment length. The compartment- 
based SSA can be described as follows. Starting with initial condition Ai(t) = ao,i, 
i = 1, 2, . . . , K, we perform six steps at time t: 

(a8) Generate two random numbers ri, r-i uniformly distributed in (0, 1). 
(b8) Compute propensity functions of reactions by ai — Ai(t)d for i = 1,2, ... ,K. 
Compute 

A"-l K 



«o = on (3.9) 



i=l i=2 



(c8) Compute the time at which the next chemical reaction takes place as t + r 
where r is given by (|2.29|) . 



STOCHASTIC REACTION-DIFFUSION PROCESSES 



19 



(d8) If r 2 < Y^f=^ a i/ a 0: then find j S {1, 2, . . . , K — 1} such that 

1 J_1 1 J 

r 2 > — > ai and r 2 < — > on. 
ao r-f a f-f 

i — 1 2—1 

Then compute the number of molecules at time t + r by 

A j (t + r)=A j (t)-l, (3.10) 

A J - +1 (t + T)=A J -+i(i) + l, (3.11) 

A i (t + r)=A i (t), fori^t^j + l. (3.12) 

(e8) If r 2 > X)^ 1 a i/ a o, then find j g {2, 3, ... , if} such that 

r2 ^(E a *+E a *) and r2< ^r(E a *+E 

u \ i=l i=2 / u \ <=1 <=2 / 

Then compute the number of molecules at time t + r by 

A i (t + T)=i4 i (t)-1, (3.13) 

A j _ 1 (t + r) = A j _i(t) + l, (3.14) 

Ai(t + r) = Ai(t), ioxi + i,i+3-\. (3.15) 

(f8) Continue with step (a8) for time t + r. 



The first term on the right hand side of (|3.9[) corresponds to reactions Ai — ► Ai + \ 
(jumps to the right) and the second term corresponds to reactions Ai — > Aj_i (jumps 
to the left). The time of the next chemical reaction is computed in the step (c8) using 
formula (|2.29|) derived previously. The decision about which reaction takes place is 
done in steps (d8)-(e8) with the help of random number r 2 . Jumps to the right are 
implemented in step (d8) and jumps to the left in step (e8). 

We want to redo the example from Section 13. li i.e. simulate 1000 molecules 
starting from position 0.4 mm in the interval [0, L) for L — 1 mm. We use K = 40. 
Since 0.4 mm is exactly a boundary between the 16th and 17th compartment, the 
initial condition is given by A 16 (0) = 500, A 17 (0) = 500 and A,;(0) = for i ^ 16, 
i ^= 17. As D = 10~ 4 mm 2 sec -1 , we have d = D/h 2 = 0.16 sec -1 . The numbers Ai(t), 
i = 1, . . . , K, at time t = 4 min, are plotted in Figure [3~37 a) as a histogram. This panel 
can be directly compared with Figure IX2T bV The computational intensity of SSA 
(a8)-(f8) can be decreased using the appropriate way to implement it in the computer. 
For example, only one chemical reaction occurs per time step. Consequently, only 
two propensity functions change and need to be updated in step (b8). Moreover, the 
formula (|3.9[) can be simplifed as follows 



K~l K K K 



ao = E] a i + E/ ai ~ ^ E/ ai ~ ai ~ aK = ^ E/ ^ ~ ai ~ aK = — ai — ax , 



1=1 i=2 i=l i=l 



where N — 1000 is the total number of molecules in the simulation (this number is 
conserved because there is no creation or degradation of the molecules in the system) . 
Hence, we need to recompute ao only when there is a change in a± or aK, i.e. whenever 
the boundary compartments were involved in the previous reaction. 
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(a) (b) 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



x [mm] x [mm] 

Fig. 3.3. Compartment-based SSA model of diffusion, (a) Numbers Aj(t), i = 1,2, . . . , K , at 
time t = 4 min obtained by SSA (a8)-(f8). We use d = D/h 2 = 0.16 sec -1 , K = 40 and initial 
condition Aig(0) = 500, Ai7(0) = 500 and Ai(0) = for i ^ 16, i ^ 17. (b) Ten realizations of the 
simulation of an individual molecule by SSA (a8)-(f8). 



SSA (a8)-(f8) does not compute the trajectories of individual molecules. However, 
we can still compute a plot comparable with Figure [3~2l a). To do that, we repeat 
the simulation with 1 molecule instead of 1000. Then, at given time t, exactly one of 
numbers Ai(t), i = 1, 2, . . . , K, is non-zero and equal to 1. This is a position of the 
molecule at time t. Ten realizations of SSA (a8)-(f8) with one molecule released at 
0.4 mm at t = are plotted in Figure l3~3l b). This panel can be directly compared 
with Figure l3~27a). 

Let p(n,t) be the joint probability that A^t) = n iy i = 1,...,K, where we 
denoted n = [m, n 2 , . . . , n K }. Let us define operators Ri,Lj : N K — > N K by 

Ri : [m, . . . ,m,n i+ i, . . . ,n K ] -> [m, ...,rii + l,n i+ i - 1, . . . ,n K ], i = l,...,K-l, 

(3.16) 

Li : [m, . . .,rii-i,ni, . . . ,n K ] [ni, . . . -1,^ + 1,.. .,n K ], i = 2,...,K. 

(3.17) 

Then the chemical master equation, which corresponds to the system of chemical 
reactions given by (13. 8|) . can be written as follows 



ar> , \ K-l . if 



d E {fa + !) p (^ n ) - ^ p ( n )} + d E {k- + x ) p ( i i n ) - n ^ p ( n )}- 

(3.18) 



The stochastic mean is defined as the vector M(t) = [Mi, M2, . . . , M^] where 



00 00 



M 4 (t) = e ^ p ( n ^ = E E • • • E n « p ( n ' *) ( 3 - 19 ) 

n ni— 0ri2— n^— 

gives the mean number of molecules in the i-th compartment, i = 1, 2, . . . , K. To 
derive an evolution equation for the stochastic mean vector M(i), we can follow the 
method from Section 12.21 - see derivation of (|2.16[) from chemical master equation 
(|2.12|) . Multiplying (|3 . 1 8|) by rij and summing over n, we obtain (leaving the details 
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to the student) a system of equations for Mj of the form 



dt 



d{M i+1 +Mi-x-2Mi), i = 2,...,K-l, (3.20) 



^±=d(M 2 -M 1 ), ®MK= d (M K - 1 -M K ). (3.21) 

System p.20j) - (|3.2ip is equivalent to a discretization of (|3.7p provided that d — D/h 2 . 
Hence, we have derived the relation between the rate constant d in (|3.8p . diffusion 
constant D and compartment length h. The solution of (|3.7p with the Dirac-like 
initial condition at x = 0.4 mm is plotted for comparison as a red solid line in Figure 
l33Ta). 

The noise is described by the variance vector V(<) = [Vi, V%, ■ ■ ■ , Vk] where 

CO 00 00 

V i (t)=J2^i-M i (t)fPM= E E E (n t -M t (t)) 2 P(n,t) (3.22) 

n ni-0 H2-0 

gives the variance of number of molecules in the i-th compartment, i = 1,2, . . . ,K. 
To derive the evolution equation for the vector V(i), we define the matrix {Vi j} by 

Vij = E ^( n ; t) ~ MiMj, for i, j = 1, 2, . . . , K. 

11 

Using (j3T2"2l) , we obtain ^ = for i = 1,2, ...,K. Multiplying (j3~T5)) by nf and 
summing over n, we obtain 

K-l 



^ E ™? p ( n ) = d E { E fa ■ + !) p (^ n ) - E ^ p ( n )} 

n j— 1 n n 

+ d E { E ^ + !) p (^ n ) - E n fa p ( n )}- ( 3 - 23 ) 



J=2 n 



Let us assume that i = 2,...,if— 1. Let us consider the term corresponding to j = i 
in the first sum on the right hand side. We get 

E + 1) p(w - E n ^ p ( n ) = E( ?ii - l ? n * p ( n ) - E n ^ p ( n ) 



= E(- 2n ? + "*) p ( n ) = - 2V > - 2M 1 + M *- 

First, we changed indices in the first sum Riii — > n and then we used definitions (|3.19p 
and (|3.22p . Similarly, the term corresponding to j = i — 1 in the first sum on the right 
hand side of (|3 . 23[1 can be rewritten as 

E nfln^ + 1) P(iVin) - E »?n<-l ^( n ) = Et 2 ™^- 1 + n ^ P ( n ) 



= 2Vi,i-i + 2M,M i _ 1 + Mi-!. 
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Other terms (corresponding to j ^ i,i — 1) in the first sum on the right hand side of 
(|3.23p are equal to zero. The second sum can be handled analogously. We obtain 

^ rfPin) = d{2V5,<_i + IM^h-i + M<_i - 2VS - 2M? + M,| 

n 

+ d{2V^, 4+1 + 2M,Af l+1 + M i+1 - 2Vi - 2M? + Mj|. (3.24) 
Using ([332]) and (|3T2H)) on the left hand side of [pT24")l . we obtain 

+ d{2M t M l+1 + 2M i M i _ 1 - 4Mf). 



2^ [ + M,+i + M,_i + 2M, ^ (3.25) 



Vi} +d{M 2 +M 1 X, (3.26) 

Fa} + rf{Mif_i + M A -}. (3.27) 

We see that the evolution equation for the variance vector V(t) depends on the mean 
M, variance V and on non-diagonal terms of the matrix Vij. To get a closed system 
of equations, we have to derive evolution equations for Vij too. This can be done 
by multiplying (|3 . 1 8|) by riiUj, summing over n and following the same arguments as 
before. We conclude this section with some consequences of (|3.20|) - (|3.21[1 and (13.25)) — 
(|3.27[) . Looking at the steady states of equations (|3.20[) - (|3.21[) . we obtain Mj = N/K, 
i = 1, 2, . . . , K, where N is the total number of diffusing molecules. Moreover, the 
variance equations imply that V% — N/K, i — 1,2, ... ,K, at the steady state. 

4. Stochastic reaction-diffusion models. In this section, we add chemical 
reactions to both models of molecular diffusion which were presented in Section [3l 
We introduce two methods for the stochastic modelling of reaction-diffusion processes. 
The first one is based on the diffusion model from Section [3T2l the second one on the 
diffusion model from Section [3. II We explain both methods using the same example. 
Namely, we consider molecules (e.g. protein) which diffuse in the domain [0, L] with 
diffusion constant D as we considered in Section[3l Moreover, we assume that protein 
molecules are degraded (in the whole domain) and produced in part of the domain, 
i.e. we consider the chemical reactions from Sections 12.11 and 12.21 in our illustrative 
reaction-diffusion model. The model has a realistic motivation which is discussed in 
more detail later in Section[5T21 I n Section ^. 3[ we present another illustrative example 
of a reaction-diffusion process incorporating the nonlinear model (|2.27p - (l2.28p . 

4.1. Compartment-based reaction-diffusion SSA. We consider molecules 
of chemical species A which are diffusing in the domain [0, L], where L = 1 mm, 
with diffusion constant D = 10~ 4 mm 2 sec _1 . Initially, there are no molecules in 
the system. Molecules are produced in the part of the domain [0, L/5] with rate 



n 

Substituting this into (|3.24[) , we get 

dVi r 

for i = 2, . . . , K — 1. Similarly, we get 



- w = 2d\v hh , 
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(a) 



(b) 




0.4 0.6 0.8 
x [mm] 




Fig. 4.1. One realization of the Gillespie SSA (a5)-(d5) for the system of chemical reactions 
\4-l\ - ^4~3\ - Gray histograms show numbers of molecules in compartments at time: (a) t = 10 min; 
(b) t = 30 min. Solution of ^4--9\ - ^4- 10\ is plotted as the red solid line. 



k p = 0.012 (im sec . This means that the probability that a molecule is created 
in the subinterval of the length 1 ^m is equal to k p dt. Consequently, the probability 
that a molecule is created somewhere in the interval [0, L/5] is equal to k p L/5dt. 
Molecules are degraded with rate k\ — 10~ 3 sec -1 according to the chemical reaction 

(EU). 

Following Section 13.21 we divide the computational domain [0, L] into K = 40 
compartments of length h = L / K = 25 fim. We denote the number of molecules of 
chemical species A in the i-th compartment [(i — l)h, ih) by A4, i = 1, . . . , K. Then 
our reaction-diffusion process is described by the system of chemical reactions 



Ai < — A 2 

d 



d d 

— A 3 ^= 



.4 



Ki 



(4.1) 



A, 



fei 



fori = 1,2,..., K, 



(4.2) 



Ai 



fori = 1, 2, . . . , K/5. 



(4.3) 



Equation (|4.ip describes diffusion and is identical to (|3.8p . In particular, the rate 
constant d is given by d — D/h 2 . Equation l|4.2p describes the degradation of A and 
is, in fact, equation (|2.1[) applied to every compartment. Equation (|4.3[) describes 
the production of A in the first K/5 compartments (e.g. in part [0, L/5] of the 
computational domain). The rate constant k 2 describes the rate of production per 
compartment. Since each compartment has length h, we have k 2 — k p h. 

The system of chemical reactions (|4.ip ~ (|4.3p is simulated using the Gillespie SSA 
(a5)-(d5). In our case, the propensity functions of reactions in (|4.ip are given as 
Ai(t)d, the propensity functions of reactions in (|4.2p are given as Ai(t)k\ and propen- 
sity functions of reactions in (|4.3p are equal to k 2 . Starting with no molecules of A in 
the system, we compute one realization of SSA (a5)-(d5) for the system of reactions 
(|4.ip - (|4.3p . We plot the numbers of molecules in compartments at two different times 
in Figure |4~T1 
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Let p(n, t) be the joint probability that Aj(t) = n,, i = 1, . . . , K, where we use 
the notation n = [m, n 2 , ■ ■ ■ , %]. Let us define operators Ri,Li : N K — > N K by 
(|3.16|) - (|3.17[) . Then the chemical master equation, which corresponds to the system 
of chemical reactions (14. ID — (14.31) . can be written as follows 



— =d 

i=l i=2 



53 { (m + !) P(R^) - rii p(n) } + d 53 {(ni + 1) P(^ n ) - P( n ) } 

t=l i=2 
if 

+ h 53 |(rij + . . . ,rii + 1, . . . .n^) - n,p(n)| 

2=1 

if/5 

53 {p{m, ■ ■ ■ ,ni - 1, ■ ■ • ,njr) -p(n)|. (4.4) 



i=l 
if/5 

i=i 



The first two sums correspond to diffusion (|4.ip . the third sum to degradation (|4.2|) 
and the fourth sum to production (|4. 3[) . The stochastic mean is defined as the vector 
M(t) = [Mi, M2, . . . , M/f] where M$ is given by (|3.19[) . This gives the mean number 
of molecules in the i-th compartment, i = 1, 2, . . . , K, at time t (averaged over many 
realizations of SSA (a5)-(d5)). To derive the evolution equation for the stochastic 
mean vector M(i), we can follow the method from Section T2.2I - see derivation of 
(|2.16|) from the chemical master equation (|2.12|) . Multiplying (|4.4|) by rij and summing 
over all rij, j = 1, . . . , K, we obtain (leaving the details to the student) a system of 
equations for Mj in the form 

^ = d{M 2 - M x ) + k 2 - hM u (4.5) 
° }l d(M i+1 + Mi_i — 2Mj) + k 2 - fciMj, i = 2,...,K/5, (4.6) 



dt 



d{M l+1 + Mi_ x - 2Af,) - feiMi, * = K/5 + 1, . . . , K - 1, (4.7) 



dt 

= d(Mjc_i - My) - fciMjf. (4.8) 

System (|4.5H — (|4.8[) is a discretized version of the reaction-diffusion equation 

da d 2 a 

-Qi= D Q^ + k mo,L/5]-kia (4.9) 

with zero-flux boundary conditions 

da , , da , , , 

Here, X[o,L/5] 1S the characteristic function of the interval [0, L/5]. Using initial con- 
dition a(-,0) = 0, we computed the solution of (|4.9[) — (|4.10|) numerically. It is plotted 
as a red solid line in Figure |4~T1 for comparison. 

The concentration of molecules in the i-th compartment is defined as Mi = Mi/h, 
i = 1, . . . , K. Dividing (|4.5|) - (|4.8|) by h, we can write a system of ODEs for Mj. It is 
a discretized version of the reaction-diffusion equation 

da d 2 a 

fa = D + k pX[o,L/5] - kia (4.11) 
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where a = a(x, t) is the concentration of molecules of A at point x and time t. The 
equation (|4. 1 1|) provides a classical deterministic description of the reaction-diffusion 
process. Its parameters D, k p and k\ are independent of h. Solving (|4.1ip is equivalent 
to solving (|4.9|) . Consequently, the red solid line in Figure [47T1 can be also viewed as a 
plot of ah where a is the solution of the classical deterministic model (|4. 1 1 [) with the 
zero- flux boundary conditions. 

4.2. Reaction-diffusion SSA based on the Smoluchowski equation. In 

this section, we present a SSA which implements the Smoluchowski model of diffusion 
from Section [3. II that is, we follow the trajectories of individual molecules. Diffusion 
of each molecule is modelled according to the model (a7)-(c7). We explain the SSA 
using the reaction-diffusion example from Section ^. II Choosing a small time step At, 
we perform the following three steps at time i: 

(a9) For each molecule, compute its position at time t + At according to steps 
(a7)-(c7). 

(b9) For each molecule, generate a random number r\ uniformly distributed in 
the interval (0, 1). If t\ < k\ At, then remove the molecule from the system. 

(c9) Generate a random number r-i uniformly distributed in the interval (0,1). 
If 7'2 < k p L/5At, then generate another random number 7'3 uniformly dis- 
tributed in the interval (0, 1) and introduce a new molecule at position r^L/b. 
Continue with step (a9) for time t + At. 

The degradation of molecules is modelled by step (b9). Equation (|2.ip implies that 
ki dt is the probability that a molecule is degraded in the time interval [t, t + dt) 
for infinitesimally small dt. SSA (a9)-(c9) replaces dt by the finite time step At 
(compare with SSA (al)— (bl)) which has to be chosen sufficiently small so that 
ki At -C 1. Similarly, the probability that a molecule is created in [0, Z/5] in time 
interval [t,t + dt) is equal to k p L/5dt. Consequently, we have to choose At so 
small that k p L/hAt is significantly less than 1. We choose At = 10 -2 sec. Then 
k\ At — 10~ 5 and k p L/5 At — 2.4 x 10~ 2 for our parameter values ki = 10~ 3 sec -1 , 
k p = 0.012 /im -1 sec -1 and L = 1 mm. Starting with no molecules of A in the sys- 
tem, we compute one realization of SSA (a9)-(c9). To visualise the results, we divide 
the interval [0, L] into 40 bins and we plot the numbers of molecules in bins at time 
10 minutes in Figure l4~2T a). The same plot at time 30 minutes is given in Figure 
I4.2f b). We used the same number of bins to visualise the results of SSA (a9)-(c9) 
as we used previously in the compartment-based model. Thus Figure 14.21 is directly 
comparable with Figure 14.11 We also plot the solution of (|4.9[) - (|4.10[) as a red solid 
line for comparison. 

4.3. Reaction-diffusion models of nonlinear chemical kinetics. In the 

previous sections, we studied an example of a reaction-diffusion model which did not 
include the second-order chemical reactions ()2.26|) . We considered only production 
and degradation, i.e. we considered chemical reactions from Sections 12.11 and 12.21 
In this section, we discuss generalisations of our approaches to models which involve 
second-order chemical reactions too. Our illustrative example is a reaction-diffusion 
process incorporating the nonlinear model (|2.2T|) — (|2.28|) . The second-order chemical 
reactions (|2.26|) require that two molecules collide (be close to each other) before 
the reaction can take place. The generalisation of SSA (a9)-(c9) to such a case is 
nontrivial and we will not present it in this paper (it can be found in pQ). Application 
of the Gillespie SSA (a5)-(d5) is more straightforward and is presented below. 
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(a) (b) 
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Fig. 4.2. One realization of SSA (a9)-(c9). Dividing domain [0,L] into 40 bins, we plot the 
number of molecules in each bin at time: (a) t = 10 min; (b) t = 30 min. Solution of $4-9\ - ^j-10\ 
is plotted as the red solid line. 



We consider that both chemical species A and B diffuse in the domain [0, L] , where 
L = 1 mm, with diffusion constant D = 10~ 4 mm 2 sec -1 . Following the method of 
Section I4TTI we divide the computational domain [0, L] into K = 40 compartments of 
length h = L /K = 25 /im. We denote the number of molecules of chemical species 
A (resp. B) in the i-th compartment \(i — l)h,ih) by Ai (resp. Bi), i = 1,...,K. 
Diffusion corresponds to two chains of "chemical reactions" : 

d d d d 

Ai^ A 2 ^ A 3 ^ . . . ^ A K (4.12) 

d d d d 



d d d d 

B 1 ^ B 2 ^ B 3 ^ ... ^ B K (4.13) 

d d d d 

Molecules of A and B are assumed to react according to chemical reactions (12.27)) in 
the whole domain with rate constants k\ = 10 _3 sec _1 and k 2 = 10 _2 sec _1 per one 
compartment, that is, 

Ai+Ai^®, A l + B l ^H), fort = 1,2,..., Jf. (4.14) 

Production of chemical species (|2.28[) is assumed to take place only in parts of the 
computational domain [0, L\. Molecules of chemical species A (resp. B) are assumed 
to be produced in subinterval [0, 9L/10] (resp. [2L/5,L]) with rate £3 = 1.2 sec -1 
(resp. fc4 = 1 sec -1 ) per one compartment of length h, that is, 

Ai, for i = 1, 2, . . . , 9JT/10, (4.15) 



B h tovi = 2K/5,...,K. (4.16) 

Starting with no molecules in the system at time t = 0, we present one realization 
of the Gillespie SSA (a5)-(d5) applied to the chemical system (|4.12j) - (|4.16p in Figure 
14.31 We plot the numbers of molecules of A and B at time 30 minutes. 



STOCHASTIC REACTION- 



DIFFUSION PROCESSES 



27 




Fig. 4.3. One realization of the Gillespie SSA (a5)-(d5) for the system of chemical reactions 
\4-12\ Numbers of molecules of chemical species A (left panel) and B (right panel) in 
compartments at time 30 minutes (gray histograms) . Solution of \4-17\ ~ \4--19\ is plotted as the red 
solid line. 



We already observed in Section 12.31 that the analysis of the master equation for 
chemical systems involving the second order reactions is not trivial. It is not possible 
to derive the equation for stochastic means as was done in Section 14.11 for the linear 
model. Hence, we will not attempt such an approach here. We also observed in 
Section fPl that the equation for the mean vector (|4.5|) — (|4.8[) was actually equal to a 
discretized version of the reaction-diffusion equation (|4.9[) - (|4.10[) which would be used 
as a traditional deterministic description. When considering the nonlinear chemical 
model (|4.14p - (|4.16p . we cannot derive the equation for the mean vector but we can 
still write a deterministic system of partial differential equations (PDEs). We simply 
add diffusion to the system of ODEs (|2~33l) - (|2~34]) to obtain 

2km 2 ~ k 2 ab + k 3 X[o,9L/io], (4.17) 
k 2 ab + hX[2L/5,L], (4.18) 
and couple it with zero-flux boundary conditions 

Iw-lw-im-iw-^ («»> 

Using initial condition a(-.O) = and &(-,0) = 0, we can compute the solution of 
(|4.17p - (|4.19[) numerically. It is plotted as a red solid line in Figure [4~3l for comparison. 
We see that (j4. 1 7[) (|4. 19[) gives a reasonable description of the system when comparing 
with one realization of SSA (a5)-(d5). However, let us note that solution of (|4.17jl - 
(I4.19| is not equal to the stochastic mean. 

The equations (|4. 1 7[) — (|4. 19|) can be also rewritten in terms of concentrations a = 
a/h and b = b/h as we did in the case of equations (|4.9[) and (|4. llf) . Let us note that 
the rate constants scale with h as k\ = kx/h, k% = k 2 /h, &3 = k^h 1 ki = kih where 
ki, fe, &3, ki are independent of h. Consequently, the equations for concentrations a 
and b are independent of h. They can be written in terms of the parameters D, k±, 
&2, ^3 and ki only (compare with (|4.1ip ). 



da da 

dt dx 2 

db _ d 2 b 

dt dx 2 
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Finally, let us discuss the choice of the compartment length h. In Sections 13.21 and 
14.11 we considered linear models and we were able to derive the equations for the mean 
vectors (e.g. (I4.5j) - (l4.8p 1. Dividing (|4z. 5[) — (|4. 8[) by h and passing to the limit h — > 0, 
we derive the corresponding deterministic reaction-diffusion PDE (|4.11[) which can be 
viewed (for linear models) as an equation for the probability distribution function of 
a single molecule (i.e. the exact description which we want to approximate by the 
compartment-based SSA). Consequently, we can increase the accuracy of the SSA by 
decreasing h. Considering the nonlinear model from this section, the continuum limit 
h — > is not well-defined. The compartment-based SSA is generally considered valid 
only for a range of h values (i.e. the length h cannot be chosen arbitrarily small); 
conditions which the length h has to satisfy are subject of current research - see e.g. 

Section 3.5]. 

5. Two important remarks. We explained SSAs for chemical reactions, molec- 
ular diffusion and reaction-diffusion processes in the previous sections. This final 
section is devoted to two important questions: 

(a) Why do we care about stochastic modelling? The answer is given in Section 
15.11 where we discuss connections between stochastic and deterministic modelling. In 
particular, we present examples where deterministic modelling fails and a stochastic 
approach is necessary. We start with a simple example of stochastic switching between 
favourable states of the system, a phenomenon which cannot be fully understood 
without stochastic modelling. Then we illustrate the fact that the stochastic model 
might have qualitatively different properties than its deterministic limit, i.e. the 
stochastic model is not just "equal" to the "noisy solution" of the corresponding 
deterministic equations. We present a simple system of chemical reactions for which 
the deterministic description converges to a steady state. On the other hand, the 
stochastic model of the same system of chemical reactions has oscillatory solutions. 
Finally, let us note that stochasticity plays important roles in biological applications, 
see e.g. [32l El [30]. 

(b) Why do we care about reaction-diffusion processes? The answer is given in 
Section I5T21 where we discuss biological pattern formation. Reaction-diffusion models 
are key components of models in developmental biology. We present stochastic ana- 
logues of two classical pattern forming models. The first one is the so-called French 
flag problem where we re-interpret the illustrative example from Sections 14. II and 14. 21 
Then we present the reaction-diffusion pattern forming model based on the so-called 
Turing instability. 

5.1. Deterministic vs. stochastic modelling. The models presented so far 
have one thing in common. One could use the deterministic description (given by 
ODEs or PDEs) and one would obtain a reasonable description of the system. In 
Sections O [2J2] O [3J2 O and [0] we studied linear models. We showed that 
the evolution equations for the stochastic mean are equal to (the discretized versions 
of) the corresponding deterministic differential equations. In Sections 12.31 and 14.31 we 
presented nonlinear models. We were not able to derive equations for the stochastic 
mean. However, we solved numerically the corresponding systems of deterministic 
equations (ODEs ([2~33l) - (j2~34]) in Section [23] and PDEs lj4l7| - l[4~T9l) in Section [43]) 
and we obtained results comparable with the SSAs, i.e. results of the SSAs looked like 
"noisy solutions" of the corresponding differential equations. Here, we discuss exam- 
ples of problems when SSAs give results which cannot be obtained by corresponding 
deterministic models. 

Let us consider the model from Section [2.31 Its deterministic description is given 
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(a) (b) 




time [min] time [min] 



Fig. 5.1. Simulation of l|5.-?j l. One realization of SSA (a5)-(d5) for the system of chemical 
reactions 1 15. j| | (blue line) and the solution of the deterministic ODE 1 15. 21 1 (red line), (a) The 
number of molecules of A as a function of time over the first two minutes of simulation, (b) Time 
evolution over 100 minutes. 



by the system of ODEs (|2.33p ~ (|2.34p . Such a system has only one nonnegative (stable) 
steady state for our parameter values, namely a s = b s = 10. It can be observed from 
Figure l2~3l that solutions of (|2.33p - (|2.34[) converge to a s and b s as t — > oo. This is 
true for any nonnegative initial condition. The results of SSAs show fluctuation about 
the means, which are roughly equal to a s and b s (they are 9.6 for A and 12.2 for B). 
However, there are chemical systems which have two or more favourable states, so that 
the corresponding ODEs have more than one nonnegative stable steady state. For 
example, let us consider the system of chemical reactions for chemical A introduced 
by Schlogl [36] 

2A ^ 3A, ^ A. (5.i) 

The corresponding ODE is given as follows 

— = - k 2 a 3 + ki a 2 - fc 4 a + fc 3 . (5.2) 
at 

We choose the rate constants as follows: fci = 0.18 min -1 , k 2 = 2.5 x 10~ 4 min -1 , 
k 3 = 2200 min -1 and fc 4 = 37.5 min -1 . Then the ODE (jO)) has two stable steady 
states a s x = 100 and a s i = 400 and one unstable steady state a u = 220. The 
solution of (|5.2[) converges to one of the steady states with the choice of the steady 
state dependent on the initial condition. Let us consider that there are initially no 
molecules of A in the system, i.e. A(0) = 0. The solution of (|5.2p is plotted in Figure 
15. If a) as a red line. We see that the solution of (|5.2p converges to the stable steady 
state a s i = 100. This is true for any initial condition A(0) £ [0,a u ). If A(0) > a u , 
then the solution of (|5.2p converges to the second stable steady state a s2 — 400. Next, 
we use the Gillespie SSA (a5)-(d5) to simulate the chemical system (|5.ip . Starting 
with no molecules of A in the system, we plot one realization of SSA (a5)-(d5) in 
Figure I5.1f a) as a blue line. We see that the time evolution of A given by SSA 
(a5)-(d5) initially (over the first 2 minutes) looks like the noisy solution of (|5.2p . 
However, we can find significant differences between the stochastic and deterministic 
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Fig. 5.2. Self-induced stochastic resonance, (a) One realization of SSA (a5)-(d5) for the 
system of chemical reactions 1 15. 311 (blue line) and solution of the deterministic ODEs (red line), (b) 
Comparison of the stochastic and deterministic trajectories in the (A, B)-plane. Nullclines of the 
deterministic ODEs are plotted as green lines. 



model if we observe both models over sufficiently large times - see Figure 15. f f b) 
where we plot the time evolution of A over the first 100 minutes. As expected, the 
solution of the deterministic model (|5.2p stays forever close to the stable steady state 
a s i = 100. The number of molecules given by the stochastic model initially fluctuates 
around one of the favourable states of the system (which is close to a s \ = 100). 
However, the fluctuations are sometimes so strong that the system spontaneously 
switches to another steady state (which is close to a S 2 = 400). This random switching 
is missed by the deterministic description. If one wants to find the mean switching 
time between favourable states of the system, then it is necessary to implement SSAs. 
Random switching between states has been found in gene regulatory networks [15U21] . 
Theoretical or computational approaches for the analysis of suitable stochastic models 
are given in [35] M ■ 

Our next example is a nonlinear system of chemical equations for which the 
stochastic model has qualitatively different behaviour than its deterministic coun- 
terpart in some parameter regimes. The presented phenomenon is sometimes called 
self-induced stochastic resonance [27] . Following an example from [5] , we consider the 
system of chemical reactions introduced by Schnakenberg [37] 

2A + B^3A, ^Z? A, B, (5.3) 

where we choose the rate constants as k\ = 4xl0 -5 sec -1 , &2 = 50 sec -1 , k$ = 
10 sec -1 and k^ = 25 sec -1 . We use the Gillespie SSA (a5)-(d5) to simulate the time 
evolution of this system. To do that, let us note that the propensity function of the 
first reaction is equal to A(t)(A(t) — l)B(t)ki. We also derive and solve the deter- 
ministic system of ODEs corresponding to (|5.3p . Using the same initial conditions 
[A, B] — [10, 10], we compare the results of the stochastic and deterministic models 
in Figure l5T27 a). We plot the time evolution of A(t). We see that the solution of the 
deterministic equations converges to a steady state while the stochastic model has os- 
cillatory solutions. Note that there is a log scale on the A-axis - numbers of A given 
by the (more precise) SSA vary between zero and ten thousand. If we use a linear 



STOCHASTIC REACTION-DIFFUSION PROCESSES 



31 
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(b) 
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Fig. 5.3. French flag problem, (a) Deterministic model, (b) Stochastic model. 



scale on the A-axis, then the low molecular fluctuations would be invisible and the 
solution of the SSAs would look as if there were "almost deterministic oscillations" , 
although it is the intrinsic noise which makes the oscillations possible. To understand 
this behaviour better, we plot the stochastic and deterministic trajectories in the 
(A, B)-plane in Figure [?^f bV We include the nullclines of the deterministic system 
of ODEs (green lines). We see that the deterministic system follows a stable nullcline 
into the steady state (red circle). The stochastic model also initially "follows" this 
nullcline (with some noise) but it is the intrinsic noise which makes it possible for the 
stochastic model to leave the stable nullcline and oscillate (again we use a log scale 
on the A-axis). 

5.2. Biological pattern formation. Reaction-diffusion processes are key el- 
ements of pattern forming mechanisms in developmental biology. The illustrative 
example from Sections l4.ll and l4.2l was a caricature of more complicated morphogen- 
esis applications [38l [33] where one assumes that some prepatterning in the domain 
exists and one wants to validate the reaction-diffusion mechanism of the next stage 
of the patterning of the embryo. In our example, we considered a chemical A which 
is produced in part [0, L/5] of domain [0, L]. Hence, the domain [0,L] was divided 
into two different regions (prepatterning) [0, L/5] and [L/5,L]. The simplest idea 
of further patterning is the so-called French flag problem [42] , We assume that the 
interval [0, L] describes a layer of cells which are sensitive to the concentration of 
chemical A. Let us assume that a cell can have three different fates (e.g. different 
genes are switched on or off) depending on the concentration of chemical A. Then the 
concentration gradient of A can help to distinguish three different regions in [0, L] - 
see Figure 15.31 If the concentration of A is high enough (above a certain threshold) , 
a cell follows the first possible program (denoted blue in Figure 15. 3p . The "white 
program" (resp. "red program") is followed for medium (resp. low) concentrations of 
A. The deterministic version of the French flag problem is presented in Figure l5~3T a). 
We consider a solution of (|4.9[) (|4. 10|) at time 30 minutes which is the red curve in 
Figure l4~TT b) or Figure [4~27 b). The solution of (|4.9|) (|4. 10[) is decreasing in space. 
Introducing two thresholds, we can clearly obtain three well-defined regions as seen in 
Figure l5~37 a). The stochastic version of the French flag problem is presented in Figure 
I5.3f b). We take the spatial histogram presented in Figure POf b). We introduce two 
thresholds (80 and 30 molecules) as before and replot the histogram using the corre- 
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(a) (b) 




x [mm] x [mm] 



Fig. 5.4. Turing patterns, (a) Numbers of molecules of chemical species A in each compartment 
at time 30 minutes; (b) the same plot for chemical species B. 



sponding colours. Clearly, the resulting "French flag" is noisy. Different realizations 
of the SSA would lead to different noisy French flags. The same is true for the SSA 
from Figure I4.1f b) . 

Our second example of patterning in developmental biology are the so-called 
Turing patterns [HJ [T71 HH1 • They do not require any prepatterning. Molecules are 
subject to the same chemical reactions in the whole domain of interest. For example, 
let us consider a system of two chemical species A and B which react according to 
the Schnakenberg system of chemical reactions (|5.3p . Let us choose the values of rate 
constants as k\ — 10 -6 sec -1 , k 2 = 1 sec" 1 , fc 3 = 0.02 sec -1 and fc 4 = 3 sec -1 . 
The corresponding deterministic system of ODEs for (|5.3|) has one nonnegative stable 
steady state equal to a s = 200 and b s — 75 molecules. Introducing diffusion to the 
model, one steady state solution of the spatial problem is the constant one (a s , b s ) 
everywhere. However, such a solution might not be stable (i.e. might not exist in 
reality) if the diffusion constants of A and B differ significantly. We choose Da = 
10 -5 mm 2 sec -1 and D B = 10 -3 mm 2 sec -1 , i.e. D b /Da = 100. To simulate the 
reaction-diffusion problem with the Schnakenberg system of chemical reactions (|5.3p , 
we follow the method of Section 14.11 We divide the computational domain [0, L] 
into K = 40 compartments of length h = L /K = 25 fim. We denote the number of 
molecules of chemical species A (resp. B) in the i-th compartment [(i— l)h, ih) by A^ 
(resp. Bi), i = 1, . . . , K. Diffusion is described by two chains of chemical reactions 
(|4.12|) - (|4.13[) where the rates of "chemical reactions" are equal to g?a = D^/Zi 2 for 
chemical species A and ds = Ds/h 2 for chemical species B. Chemical reactions 
(|5.3|) are considered in every compartment (the values of rate constants in (|5.3|) are 
already assumed to be expressed in units per compartment). Starting with a uniform 
distribution of chemicals Aj(0) = a s = 200 and Bi(0) — b s — 75, i — 1,2, . . . , K, 
at time t = 0, we plot the numbers of molecules in each compartment at time t = 
30 minutes computed by SSA (a5)-(d5) in Figure 15.41 To demonstrate the idea of 
patterning, compartments with many molecules (above steady state values a s or b s ) 
are plotted as blue; other compartments are plotted as red. We see in Figure HT47 a) 
that chemical A can be clearly used to divide our computational domain into several 
regions. There are two and half blue peaks in this figure. The number of blue peaks 
depends on the size of the computational domain [0, L] and it is not a unique number 
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in general. The reaction-diffusion system has several favourable states with a different 
number of blue peaks. As discussed in Section I5TT1 the solution of the corresponding 
deterministic model converges to one of the favourable (stable steady) states of the 
system. The stochastic model enables stochastic switching between the favourable 
states, i.e. between the states with a different number of blue peaks. 

6. Discussion. We presented SSAs for systems of chemical reactions and molec- 
ular diffusion. Then we presented methods for simulating both reactions and diffusion 
at the same time. The algorithms for simulating (spatially homogeneous) systems of 
chemical reactions were based on the work of Gillespie [IB] . We did not focus on the 
computer implementation of the algorithms. We chose simple examples which can 
be simulated quickly. If one considers systems of many equations, there are ways 
to make the Gillespie SSA more efficient [16]. For example, it would be a waste of 
time to recompute all the propensity functions at each time step. We simulate one 
reaction per one time step. Therefore, it makes sense to update only those propensity 
functions which are changed by the chemical reaction which was selected in step (d5) 
of SSA (a5)-(d5). 

We only briefly touched on the concept of the Fokker-Planck equation [33] when 
we discussed the Smoluchowski description of diffusion. It is worth noting that there 
are interesting connections between the chemical master equation (which is equivalent 
to the Gillespie SSA) and the Fokker-Planck equation which gives the time evolution 
of the probability distribution. Such connections are discussed (through the so-called 
chemical Langevin equation) in [20 . The Smoluchowski equation is actually the same 
mathematical object as the chemical Langevin equation, i.e. the stochastic differential 
equation [2]. An algorithmic introduction to stochastic differential equations can be 
found in [25] . 

We presented two models of diffusion in this paper. One was based on the chain of 
"chemical reactions" Q3.8P computing the time evolution of the numbers of molecules 
in compartments. Coupling this model with the modelling of chemical reactions is 
straightforward and presented in Section 14. l\ such a compartment-based approach 
is used e.g. in [40] [231 [22] • The second model for molecular diffusion was based on 
the Smoluchowski equation (|3.6[) . It was an example of the so-called position jump 
process, that is, a molecule jumps to a different location at each time step. As a result, 
the trajectory of a molecule is discontinuous. The individual trajectories of diffusing 
molecules can be also modelled using the so-called velocity jump processes [29] . that 
is, the position of a molecule x(t) follows the deterministic equation dx/dt = v where 
v(t), the velocity of the molecule, changes stochastically. Such stochastic processes can 
be used not only for the simulation of diffusing molecules but also for the description 
of movement of unicellular organisms like bacteria (TOl [11] or amoeboid cells j!2j . 
Velocity jump processes can be also described in terms of PDEs for the time evolution 
of the probability distributions to find a particle (molecule or cell) at a given place. 
Such equations are not exactly equal to the diffusion equation. However, they can be 
reduced in the appropriate limit to the diffusion equation [431 [7] . A classical review 
paper on diffusion and other stochastic processes was written by Chandrasekhar [3], 
a nice introduction to random walks in biology is the book by Berg [3J . 

In this paper, we used only reflective boundary conditions, that is, particles hitting 
the boundary were reflected back. Such boundary conditions are suitable whenever 
there is no chemical interaction between molecules in the solution and the boundary 
of the domain. Considering biological applications, it is often the case that molecules 
(e.g. proteins) react with the boundary (e.g. with receptors in the cellular mem- 
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brane). Then the boundary conditions have to be modified accordingly. It has to be 
assumed that some molecules which hit the boundary are reflected and some molecules 
are adsorbed by the boundary (e.g. become bound to the receptor or take part in 
membrane-based chemical reactions). The probability that a molecule is adsorbed 
rather than reflected depends on the chemical properties of the boundary and also on 
the SSA which is used for modelling (further details are given in [318]). 

Our analysis of SSAs was based on the chemical master equation. We successfully 
derived equations for the means and variances in illustrative examples which did 
not include second-order reactions. Other first-order reaction networks can be also 
analysed using this framework |13j . The nonlinear chemical kinetics complicates the 
mathematical analysis significantly. We can write a deterministic description but it 
might be too far from the correct description of the system [35] . A review of more 
computational approaches for the analysis of SSAs can be found in [26) . Applications 
of such methods to stochastic reaction-diffusion processes is presented in [31] . 
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